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With the help of the improved Monte Carlo renormalization-group 
scheme, we numerically investigate the renormalization group flow of the anti- 
ferromagnetic Heisenberg and XY spin model on the stacked triangular lattice 
(STA-model) and its effective Hamiltonian, 2A^-component chiral (f)^ model 
which is used in the field-theoretical studies. We find that the XY-STA 
model with lattice size 126 x 144 x 126 exhibits clear first-order behavior. 
We also find that the renormalization-group flow of the STA model is well 
reproduced by the chiral cf)^ model, and that there are no chiral fixed points 
of renormalization-group flow for N = 2 and 3 cases. This result indicates 
that the Heisenberg- STA model also undergoes first-order transition. 

PACS numbers: 75.10.Nr 



I. INTRODUCTION 

The critical behavior of the antiferromagnetic vector spin models on the stacked triangu- 
lar lattice (STA) is still a controversial issue even after twenty years of extensive studies by 
means of experimental, field-theoretical, and numerical methods. See0'ii'@ for recent works, 
and! for a review. 

The field-theoretical renormalization group (RG) analysis tells that when the number 
of spin component is greater than some threshold value N^-, there are so-called "chiral 
fixed points" of the RG which control the critical behavior and are characterized by novel 
values of critical exponents, while for N < Nc such fixed points disappear and the phase 
transition is of first order (seeifor a review). Theoretical estimations of Nc have been made 
by various authors by means of e-expansionii'Q, fixed dimensional perturbationiJl, local 
potential approximation!, and effective average action approach^, but the estimated values 
range from negative value to 6.5, depending on the method employed (however, recent results 
tend to be around 6.0). Thus the critical behaviors of the physical relevant cases N = 2 (XY 
spin) and = 3 (Heisenberg spin) are still unclear. A number of experimental studies! and 
numerical simulation^llll have yielded results which suggest second-order phase transition 
for N = 2 and 3. However, recently it has been pointed out that the critical exponent 
7], calculated from the scaling relation, becomes negative in some of these results, which is 
unphysical and indicates that the observed critical behavior is in fact pseudo-critical behavior 
induced by the slow RG fiowEl&i. 

The present work intends to clarify the issue by numerically observing the RG flow and 
investigating whether there are chiral fixed points or not for several values of A^. The paper is 
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organized as follows: in the next section, we will present the method by which we observe the 
RG flow in the Monte Carlo simulations. Details of the Monte Carlo simulations are given 
in Sec. |ITT[ In Sec. the RG flow diagrams obtained from the simulations for iV = 2, 3, 8 
cases are presented. The last section is devoted to the concluding remarks. 



II. NUMERICAL OBSERVATION OF THE RG FLOW 
A. chiral 0^ model 

The critical behavior of the A^-component STA model is essentially described by the 
following 2A^-component Ginzburg-Landau- Wilson Hamiltoniani: 



H = Jdx\K + +r{$l + $1) + u ($1 



where = (0a, 0b) is + A^-component vector field defined on the continuum space. In the 
field-theoretical studies, K is fixed to unity to eliminate ambiguity of coefficients induced 
by a trivial rescaling of 0, namely, — >■ c0. In this regularization scheme, r plays a role of 
temperature. In the present work we concentrate on the case f > 0. Figure |l| depicts the 
renormalization group flow of the scaling variables r, m, and v. There are some trivial fixed 
points, namely: 

H : high-temperature fixed point (r, M,t>) = (oo,0,0), 

LO : anisotropic low-temperature fixed point (—00,00,00), 

LI : isotropic low-temperature fixed point (— 00, oo,0), 

G : Gaussian fixed point (0,0,0). 

The flow on the critical plane projected onto the u-v plane is shown in Fig. |^ for (a) 
N < Nc and (b) N > Nc cases. When N > Nc, there are stable and unstable chiral fixed 
points denoted by Ci and C2, respectively, beside the 0{2N) symmetric and the Gaussian 
fixed points denoted by O and G, respectively. The two fixed points Ci and C2 approach as 
decreases, and annihilate each other at A^ = N^. 

In the isotropic 0^ models, the correction-to-scaling term (distance to the Wilson-Fisher 
fixed point) rapidly converges to zero as where 6 is a renormalization factor and uj ~ 0.8 
is the correction-to-scaling exponent. Therefore the asymptotic critical behavior can be 
observed in finite systems with moderate size. However, when an anisotropic quartic term 
is included in the Hamiltonian, it in general induces very slow RG flow along the anisotropy 
direction. The correction-to-scaling exponent at an anisotropic fixed point (if it exists) is 
usually of order 0.1 or less, and extremely large system is needed to observe asymptotic 
critical behavior in the simulation of finite systems. Even when there are no chiral fixed 
points and the RG flow eventually diverges, the RG flow in the intermediate region is very 
slowl and it is very hard to observe asymptotic first-order behavior. Thus, in the simulation 
of finite systems, it is crucial to check the convergence to the fixed point of the renormalized 
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quartic coupling constants u and v. Now let us consider how to observe these renormalized 
coupling constants in the numerical simulations. In the continuum theory, they describe the 
behavior of renormalized field variables 0(x; /) defined as follows: 

0(x; l)= [ exp(ik • x)0(k)rfk, (2) 

"'|k|<l/i 

— # — * 

where / is a cut-off length and 0(k) denotes Fourier component of the field 0. We denote by 
Ui and Vi the renormalized coupling constants of cut-off length /: the renormalization flow 
in Fig. ^ is a trajectory of ui and vi on the critical plane when / is increased. The renor- 
malization behavior of these quantities may be well observed via their conjugate quantities, 
which are easier to observe in numerical simulations: 

((0'(x;/))^), (3) 



(0^(x; 1) - (0,(x; /) • 0,(x; l)f). (4) 

Of course, these quantities are affected not only by ui and vi but also by other irrelevant 
scaling flelds, such as coefficients of higher order terms like |0|^ and However, these 
irrelevant scaling fields rapidly vanishes (faster than the leading correction-to-scaling term in 
the isotropic model, ~ Z""'^) and their effects can be ignored after sufficient renormalization 
as long as ui and vi converge (or diverge) very slowly. 

Note that the definition of the renormalized field 0(x; /) needs scaling prefactor so that 
the coefficients such as ui and vi remain finite when / goes to infinity. The following scaling 
scheme is simple and suitable for numerical simulationsEl: 

b (x; /) = (5) 

'(0(x;/)2) 



Substituting 0(x; /) by 0'(x; /), equation (|]) and (^) lead as follows: 

((0^(x;/))^) 

((02(x;O))^' 



(6) 



(0g(x;O0g(x;O-(0,(x;O-0,(x;/))^) 

Basically, in the present work, we observe these quantities in the numerical simulations 
of finite lattice version of the Hamiltonian (|l|) defined as follows: 
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(8) 



where = (J)a{i), 0b(«)) is an + A^-component vector defined on the lattice site i of an 
L X L X L cubic lattice with periodic boundary condition being imposed, the summation 
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J2<ij> runs over all nearest neighbor pairs of the lattice site. We use a restriction u + v/i — 
2r = so that the minimum of the Hamiltonian takes place at \4>{i)\ = 1 to eliminate the 
ambiguity of the trivial rescaling of 0. Actually we use the following parameterization: 

r = -2A, u = X{l + A), v = A\A (9) 

where A controls the strength of the anisotropy and A controls the "hardness" of the spin — 
the larger A is, the smaller the fluctuation of \4>{i) \ is. A line A = corresponds to u—v/4 = 0, 
on which instability of the Hamiltonian (||) occurs; thus A should be positive. 

In the lattice models, renormalized field corresponds to block spin variables. For example, 
consider a block spin variable (j)b,L defined on a block of b x b x b spins in an L x L x L 
lattice model. One can calculate the quantities and (|^ using 0^^^, which we denote by 
Ub^L and Vb^L, respectively, and check the convergence or divergence of the renormalization 
flow. In the conventional MCRG scheme", one calculates and (|^ for fixed value of 
L and increasing value of b. Note that, after the "renormalization" by the block spin 
transformation, the number of spins, {L/bY, decreases and the Hamiltonian does not stay 
in the same Hamiltonian space, unlike the continuum case. Therefore the number of block 
spins must be sufficiently large, i.e. L/b ^ 1, so that (||) and (|^) take the same value at the 
fixed point for different values of b. 

The present author proposed another scheme0, in which the ratio L/b is fixed and both 
L and b are increased. Let us denote a renormalized Hamiltonian which contains block 
spins by Hi{K, r, u, v, {wi}), where K, r, u, v are the renormalized coefficients of terms which 
correspond to that in the original Hamiltonian (^, and {wi} denote set of coefficients of 
higher order irrelevant terms such as which is not included in the original Hamilto- 
nian. Now consider the behavior of (pi^i and 4'2L,2L'- the behavior of both quantities are 
described by a renormalized Hamiltonian which contains only 1 spin, but with different 
coefficients. Here, 02L,2L is a renormalized spin, obtained by applying real-space renormal- 
ization of factor-2L to the original spins. Let us regard this factor-2L renormalization as 
two successive renormalizations, namely factor-2 renormalization at first and then factor-L 
renormalization. Accordingly, the original Hamiltonian H2L{K,r,u,v, {0}) is at first renor- 
malized to Hl{K' ,r' ,u' ,v' , {wi}), then again renormalized to Hi{K" , r" , u" , v" , {w'^}) . If 
the starting Hamiltonian is at a fixed point and L is sufficiently large, coefficients before 
and after the factor-2 renormalization should remain unchanged except for the higher or- 
der coefficients At this stage, U2l,2L and V2l,2L can be calculated (in other words, 
factor-L renormalization) from Hl{K' ,r' ,u' ,v' , {wi}), while Ul,l and Vl,l are calculated 
from HL{K,r,u,v, {0}). Since the coefficients of these two Hamiltonians are the same (ex- 
cept the irrelevant ones), U2l,2L = Ul,l + 0{L^'^) and V2l,2L = Vl,l + 0{L^'^) should be 
satisfied, where u denotes correction-to-scaling exponent of the higher order terms such as 
The so-called "phenomenological renormalization group method"lli uses this property 
to locate the critical point. 

When the Hamiltonian is not at a fixed point, the coefficients before and after the first 
factor-2 renormalization differ: they moves along the renormalization ffow. This difference is 
reflected by Ul^l and If the renormalization flow slowly converges to some flxed point, 
Ul^l and Vl^l slowly converge to some flxed value, while when the phase transition is of flrst- 
order and the renormalization flow diverges, so do Ul,l and Vl.l^^- We will henceforth denote 
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these quantities simply by Ul and Vl, and investigate their behavior at the critical point 
when L is increased. Here their explicit form is written down for the readers' convenience: 

Ul = ^ ^ ^ , (10) 



' ^''^ 

where M = {MayMi,) = Yl,i'P{i)- In the past numerical studies in which only a specific 
model was investigated, the only way to reach the final fixed point (if it exists) or to ob- 
serve asymptotic first-order behavior was to simulate extremely large systems, and it was 
impossible in general. In the present work, we investigate various Hamiltonians and scan 
the Hamiltonian space in search of fixed points. This method allows us to determine the 
asymptotic critical behavior by simulations of moderately large systems. 

Note that, in the theoretical works one can concentrate on the RG flow on the critical 
plane on which inverse susceptibility vanishes, while in the numerical works one must deter- 
mine the critical point from the numerical data. For this purpose we use a quantity which 
corresponds to the following quantity in the continuum theory: 

(/|.|<iMk/OV^^(k)rf^k) ^ ^, ((V0(x;/))2) 
((02(x;/))) ((0^(x;/))) ■ 

In the lattice models, the above quantity corresponds to the normalized correlation be- 
tween two adjacent block spins. Actually we observe the following quantity: 

^ Wk.)./(-k.)> 

— * — # 

where 0(ki) = X^r 0('^) exp(iki ■ r) with ki = (27r/L, 0, 0). 

Figure |^ depicts the RG flow in Ul-Vl-Cl space: a number of trivial fixed points, 
namely, high-temperature fixed point (1 + ^'^^^^j^^i)^'' , 1), isotropic low-temperature 

fixed point (1, ^^^,0), and anisotropic low-temperature fixed point (1,|,0) are denoted 
by H, Lq, and Li, respectively. Beside these points, there are some trivial manifolds: when 
the probability distribution of the order parameter M is isotropic, a relation Vl = 4^?q^ Ul 

is satisfied, while in the strong anisotropy limit where Ma ■ = and \Ma\ = \Mf,\ hold, 
Vl = Ul/4: is satisfied. Atthe Gaussian fixed point, the behavior of the finite system is 
governed by the zero-modelij, therefore the Gaussian fixed point lies on the Cl = plane. 
The RG flow on the critical plane, projected onto the Ul-Vl plane, is shown in Fig. ^ for 
(a) N < Nc and (b) N > cases. In both cases, the RG flow along the (approximately) 
horizontal direction rapidly converges, while along the (approximately) vertical direction the 
flow is expected to slowly converge /diverge. 

To obtain the RG flow on the critical plane, we investigate the Hamiltonian (|^) for 
various values of A and A, tuning K so that Cl = C^l is satisfied and observe the difference 
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between {Ul,Vl) and {U2l,V2l)- It should be noted that the above definition of critical 
plane induces systematic deviation. For example, consider the RG fiow of the isotropic 
0"^ model (set w = in (H) ). The fiow of (Ul,Cl) near the Wilson-Fisher fixed point is 
depicted in Fig. || (a)0. If one define the critical point by Cl{K) = C2l{K), the fiow of 
Ul becomes the one shown in Fig. ^(b), thus the arrow tends to "overshoot" the RG fixed 
point. However, the direction of the fiow can be correctly estimated, and the systematic 
error vanishes as one approaches the fixed point. 



B. STA models 

We also observe the RG fiow of the following STA model: 



H = -KY. J,,S, ■ 4-, (14) 



<ij> 

where Si denotes a two-component (XY) or three-component (Heisenberg) vector spin with 
\Si\ = 1 defined on the lattice site i of the stacked triangular lattice. We set Jij = — 1 
(antiferromagnetic) for intra-plane nearest neighbor pairs, Jij = 3/4 (ferromagnetic) for 
inter-plane nearest neighbor pairs, and = otherwise, so that the Fourier transform of 
Jij near the antiferromagnetic modeS kAF = (2vr/3,0,0) becomes isotropic in the fc-space, 
i.e. J{kAF + A;) ~ Jik-AF) + c\k\'^ + 0(|/c|^). Note that the critical values of quantities such 
as f/i, Vl, and Cl do not depend on the microscopic lattice structure, but do depend on the 
macroscopic lattice structure such as boundary conditions and aspect ratios§@. Therefore 
we use rectangular system which contains x Ly x spins (see Fig. imposing periodic 
boundary condition for all three directions. The aspect ratio then becomes Lj. : : 
L^. We simulate the system with {L^,Ly,L^) = (21,24,21), (42,48,42), (84,96,84), and 
(126, 144, 126), all of which give aspect ratio 1 : 0.99 ■ ■ ■ : 1. 
The order parameter is defined as followsi: 

Ma = MS{kAF)] = Ma- \mb - \mc, (15) 

M, = lm[S{kAF)] = ^Mb - ^Mc, (16) 

where S{kAF) denotes Fourier component of Si at kAF, and Ma, Mb, Mc denote the mag- 
netization on the three sublattices. Then Ul and Vl are calculated from Eq. ( |T0| ) and Eq. 
(p!T|). Cl is calculated as follows: 

^ _ {\S{kAF + h)\' + \S{kAF-h)\') 

where ki is the smallest, non-zero momentum in the fc-space. We observed the values of 
Cl for three kinds of direction, namely, ki = (|^,0,0), (0, , 0), and (0,0, |^), and 
confirmed that they coincide each other within statistical errors. This indicates that the 
simulated system is spatially isotropic. 
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III. MONTE CARLO SIMULATION 



A. chiral cf)^ models 

We used the usual single-spin update Metropolis algorithm, since there is no spin reflec- 
tion axis which preserves the anisotropy term in (H) and the cluster algorithmiil cannot be 
applied. In a single-spin update process, the new spin value was chosen as follows: 



Mit) + -j^G (18) 




where G is an 2A^-component vector whose components are independent Gaussian random 
variables with average and variance 1. The optimum value of the amplitude of random 
move, Rg, was determined by simulating small systems. The optimum value was found to 
be Rg ~ 0.3 for all A^, A, and A. Each Metropolis sweep was followed by one overrelaxation- 
type update^. Hereafter we refer one Metropolis sweep plus one overrelaxation-type sweep 
as "MCS". 

Physical quantities were observed for every Cint x L"^ MCS, where Cint is a constant 
which is adjusted so that the correlations between successively observed quantities become 
less than 0.70, namely: 

where Xt denotes observed value of a physical quantity X at the tth observation. The value 
of CjNT varies from 0.1 to 2.0, depending mainly on A. Actual values will be summarized 
in the later sections. For each values of L, N, A, and A, physical quantities were observed 
10^ to 10^ times, depending on the accuracy required to observe the RG flow. Statistical 
errors were estimated by the Jackknife method. Histogram method^ was used to calculate 
thermal averages at K slightly away from that where the simulation is actually carried out. 

Simulations were mainly carried out on Fujitsu VPP5000 vector processors at JAERI. 
The checkerboard-wise decomposition of the lattice was used for vectorization. Furthermore, 
we simultaneously simulated a number of independent systems, whose spins were stored into 
one long vector. This promoted the vectorization, especially when L is small, and accelerated 
the calculations. For = 8 and L = 16 system, one MCS took 0.84 ms on the VPP5000. 



B. STA models 

In the frustrated spin models, the cluster update algorithm tends to flip all the spins 
and does not accelerate the simulation. Therefore we again used only the single-spin update 
Metropolis algorithm, followed by an overrelaxation sweep, i.e. 180 degree rotation of a spin 
Si with respect to its local field J2j JijSj- In the Metropolis update, the new spin direction 
was chosen as follows: 

= Vl - cos{e), Sy = Vl - sm{e), = z, (20) 

for Heisenberg case and 
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S, = cosi9), Sy = sm{9), (21) 

for XY case, where z and 6 are uniformly distributed random numbers in the range [—1,1) 
and [0,27r), respectively. Physical quantities were observed for every x Ly /20 MCS for 
the Heisenberg case and x Ly/ 10 MCS for the XY case. This interval was long enough 
to satisfy the condition (|1^) as long as the energy histogram has no double-peak. Physical 
quantities were observed 10'^ times for all system sizes up to (L^, Ly, L^) = (42,48,42). 

To observe the energy histogram of the larger system, we also simulated very large sys- 
tems (Lj;, Ly, Lz) = (84, 96, 84) and (126, 144, 126) for the XY case. Owing to the prominent 
double-peak of the energy histogram, autocorrelation time was much longer than x Ly/ 10 
MCS in these sizes. We used 4 x 10^ MCS for the (84, 96, 84) case and 7 x 10^ MCS for 
the (126, 144, 126) case, which were 1000 times longer than the auto-correlation time. One 
MCS oiN = 3 and (L^, Ly, L,) = (84, 96, 84) system took 27 ms on the VPP5000. 



C. Stiefel's Vn,2 models 

Stiefel's Vn2 nio dei corresponds to the Hamiltonian (|^) with X = oo and A = cxd, in 
which restrictions (pa - 4>b = 0, 0^ = 0^ = 1/2 are imposed to all the spins. We also simulated 
this model and investigated the RG flow. In the Monte Carlo simulation, only the single- 
spin Metropolis update was used. We only simulated the ¥3^2 niodel, since the ¥2,2 niodel is 
known to exhibit strong first-order behaviori. In the Metropolis update, the new spin value 
was chosen as follows: 

0^ = {z,Vl- z'^cos9i,Vl- z^sm9i)/V2, (22) 
(pb = (Vl — 2^ sin 6*2, — sin 61 cos 6^2 — z cos 61 sin 62, cos 6*1 cos^^2 — -2 sin 6^1 811162) /V2 (23) 

where z and 61^2 are uniformly distributed random numbers in the range [—1, 1) and [0, 2tt), 
respectively. Physical quantities were observed for every L^ Metropolis sweeps. This interval 
was long enough to satisfy the condition (|19]) for all system sizes up to L = 48. Physical 
quantities were observed 10*^ times for these system sizes. 

To observe the energy histogram of the larger system, we also simulated a very large 
system L = 64 and L = 80. Like the STA-XY case, a prominent double-peak of the energy 
histogram emerges and autocorrelation time gets much longer than L^ MCS in these sizes. 
We used 6 x 10^ MCS for the L = 64 case and 2 x 10^ MCS for the L = 80 case, which 
were about 1000 times longer than the auto-correlation time. For the largest size L = 80, 
one Metropolis sweep took 22 ms on the VPP5000. 



IV. RESULTS 
A. XY case 

Parameters A, A, and K at which simulations were carried out are summarized in Table 
|. Figure shows the histogram of energy per spin, E = —Y^ JijSiSj / L^LyL^ of the STA- 
XY model with {L^,Ly,L,) = (126,144,126) and (84,96,84) at K = 0.77262. It can be 
seen that the valley between the two peaks deepens as the system size increases. This is a 
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clear evidence of first-order transition. This peak does not appear when the size is smaller, 
therefore the double peak has not been observed in the past studies of STA-XY model in 
which smaller systems have been used. Note that the first-order behavior has already been 
observed in the Monte Carlo simulations of other models which possess the same symmetry 
as the STA-XY model, such as Stiefel model and restricted STA modefl, and quasi-one 
dimensional STA mo delli. These models have stronger anisotropy than the STA-XY model 
and the first-order behavior can be observed in smaller systems. 

The two peaks of the histogram correspond to ordered and disordered phases, and the 
energy difference between the peaks provides lower bound of the latent heat per spin A. 
On the other hand, experimental measurement of specific heat of the STA-XY materials 
such as CsMnBrs or Holmiumil indicates second-order transition or very weak first-order 
transition whose latent heat is beneath the resolution of measurements, thus the upper 
bound of the latent heat is given. To compare the strength of the first-order transition, an 
adimensional parameter c = A/fc^T^ is usually used, where ks denotes Boltzmann constant 
and Tc denotes the transition temperature. From Fig. 0, the upper bound of the latent heat 
per spin is estimated as A > 7.8 x 10~^ and c > 6.0 x 10"^, while Ref.ii provides estimates 
c < 2.9 X 10~^ for CsMnBrs and c < 2.8 x 10~^ for Holmium. These upper bound of c 
obtained from experiments is far smaller than the lower bound of c obtained in the present 
work, thus there is an apparent contradiction between them. 

Note that the value of c may depends on the strength of the easy-plane anisotropy of 
the spin and the ratio between inter and intra layer coupling constants. The result of Ref.il 
indicates that when the inter-layer coupling is much larger than the intra-layer coupling, the 
first-order transition becomes stronger. On the other hand, a weak easy-plane anisotropy 
will lead to Heisenberg-to-XY crossover which weaken the transition, since the transition 
is weakened as the number of spin component is increased and eventually becomes second- 
order. Actually, we will see in the following section that the STA-Heisenberg model does not 
exhibit as strong a transition as the STA-XY model. Thus the most plausible explanation 
of the apparent discrepancy between the experimental results and the present work may 
be that the finiteness of the easy-plane anisotropy in the experimental materials induces 
Heisenberg-to-XY crossover and the transition is weakened. 

Figure | and ^ show the KG flow of {Ul, Vl) at the strong and weak anisotropy region, 
respectively, of the N = 2 chiral 0^ model. Figure ^ also shows the same plot for the STA- 
XY model. One can see that the STA-XY model is on the "runaway" trajectory, which 
is another evidence of first-order transition. It can also be seen that there are no stable 
fixed points other than the 0(4) symmetric one. This means that any models or materials 
which possess the same symmetry as STA-XY model, namely x 0(2), exhibit first-order 
transition. 

Fig. 1^ shows the vertical velocity of the RG flow V2L — Vl plotted against Vl- The 
plots for the two different sizes L = 6 and L = 8 coincide within statistical errors. This 
indicates that the system size is large enough to eliminate finite-size artifact of the RG fiow. 
Thus we can safely extrapolate the asymptotic critical behavior from the finite-size RG fiow. 
Moreover, the plots of the chiral 0^ model and the STA-XY model are essentially the same. 
This means that the critical behavior of the STA-XY model is well reproduced by the chiral 
0^ model. 
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B. Heisenberg case 



Parameters A, A, and K at which simulations were carried out are summarized in Table 
|I[ Unlike the STA-XY case, the STA-Heisenberg model does not show the double-peak 
behavior within the simulated size {L^, Ly, L^) = (84,96,84). However, the ¥3^2 model with 



L = 80 shows a clear double-peak behavior. Figure |TT] shows the histogram of energy per 
spin, E = 3-L~^ J2 Jij{Si — SjY at K = 0.217685. Although it is necessary to simulate larger 
system and show that the valley between the two peaks deepens to confirm the first-order 
transition, this double-peak strongly indicate that the transition is of first-order. 

Figure O and O show the RG flow of {U Vl) at the strong and weak anisotropy region. 



respectively, of the = 3 chiral 0^ model. Figure ^ also shows the same plot for the STA- 
Heisenberg model and the V3,2 model. It can be seen that the STA-Heisenberg model and the 
V3,2 model are on the "runaway" trajectory, therefore undergo first-order transition in the 
thermodynamic limit. It can also be seen that there are no fixed points other than the 0(6) 
symmetric one. This means that any models or materials which possess the same symmetry 
as STA-Heisenberg model, namely Z3 x 0(3), exhibit first-order transition. From Fig. ^ we 
roughly estimate that L = 96 STA-Heisenberg model needs factor-2 renormalizations three 
or more times to reach the strong first-order region, which means that the clear first-order 
behavior will be observed in an L ~ 800 system. 



Fig. |1J shows the vertical velocity of the RG flow V2L — Vl plotted against Vl- In this 
case, this quantity depends not only on Vl but also on Ul, although the dependence is 
not so strong. Therefore the plot can not be fitted to a single line well. However, there 
is no systematic deviation between the plots of the two different sizes. Thus we can safely 
extrapolate the asymptotic critical behavior from the finite-size RG flow. As in the XY case, 
the plots of the chiral (j)^ model and the STA-Heisenberg model are essentially the same and 
the critical behavior of the STA-Heisenberg model is well reproduced by the chiral 0^ model. 

C. = 8 case 

It is necessary to demonstrate that the MCRG scheme presented in Sec. |I| is able to 
detect the chiral fixed point if it does exists. For this purpose, we investigate N = 8 case, 
which is greater than any theoretical estimate of Nc in literaturJfflM. Parameters X, A, 
and K at which simulations were carried out are summarized in Table |TT[ Figure 15 shows 



the RG flow of {Ul, Vl) for L = 8 16. The stable chiral flxed point Ci and the unstable 
chiral flxed point O2 are clearly seen. 

Fig. |16| shows the vertical velocity of the RG flow V2L — Vl plotted against Vl, using 
only the data near the stable flxed point Ci. Unlike the N = 2 and the A^ = 3 case, this 
plot has a zero. The plots of the two different size coincide within statistical error, which 
indicates that the simulated system is large enough to observe the RG flow. 

V. CONCLUSION 

The large scale Monte Carlo simulation of the STA-XY model has revealed that the 
transition is of flrst-order. This is the flrst numerical work to conflrm the flrst-order transi- 
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tion of the simple and isotropic STA-XY model. The result also indicate that the transition 
is strong enough to be observed in experiments. But the finiteness of the easy-plane spin 
anisotropy of the actual materials will weaken the transition compared to the STA-XY model 
simulated in the present work, in which the easy-plane anisotropy is infinitely strong. 

We have also investigated the RG flow the 2A^-component chiral 0^ model numerically 
and have found that there are no chiral fixed points in = 2 and A^ = 3 cases, and we 
conclude that the STA-Heisenberg model also undergoes first-order transition. By a rough 
estimate, we expect that the STA-Heisenberg model exhibits strong first-order behavior 
when the size is larger than L = 800. For N = 8 case, we have found the chiral fixed points, 
thus the value of Nc is estimated as 3 < Ac < 8. Further study to determine the precise 
value of Nc may be of some interest, because such a study will serve to test the accuracy of 
field theoretical techniques. However, once the nature of physical relevant case A^ = 3 and 
A^ = 2 is clarified, the precise value of Nc is of less importance. 
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H v=0 v=oo 




FIG. 1. Renormalization flow of r, u, and v. H, G, Lq, and Li denote the high-temperature, 
Gaussian, isotropic low-temperature, and anisotropic low-temperature fixed point, respectively. 



(a) N<Nc 



(b) N>Nc 

>\7 



FIG. 2. Renormalization flow on the critical plane projected onto the u-v plane for (a) < A^c 
and (b) N > Nc- Shaded area is unstable region. G, O, Ci, and C2 denote Gaussian, 
0(2A?^)-isotropic, stable chiral, and unstable chiral fixed points, respectively. 
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L 4 (N+1) L 



UJkJ 



FIG. 4. Renormalization flow on the critical plane projected onto the Ul-Vl plane for (a) 
N < Nc and (b) N > Nc cases. The symbol O, Ci, and C2 denote 0(2A?^)-isotropic, stable chiral, 
and unstable chiral fixed point, respectively. 




FIG. 5. (a) Renormalization flow of Cl and Ul for isotropic model v = near the Wilson-Fisher 
flxed point, (b) RG flow of Ul{Kc) tends to "overshoot" the fixed point (see the main text). 



LX 




FIG. 6. Rectangular system used for the simulation of stacked triangular antiferromagnet. The 
figure shows a system with Lx = 6,Ly = 8. This layer is stacked Lz times. 
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£::Energy per spin 

FIG. 7. Histogram of en- 

ergy per spin of the STA-XY model with the size {L^,Ly,L^) = (84,96,84) and (126,144,126) 
&%K = 0.77262. 
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1.1 1.15 1.2 1.25 1.3 1.35 1.4 1.45 



FIG. 8. RG flow of tlie N = 2 cliiral c/)'^ model and the STA-XY model at the strong anisotropy 
region. The numbers in parenthesis are the values of (A, A) of the chiral model. 
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FIG. 9. RG flow of the N = 2 chiral (f>^ model at the weak anisotropy region. The numbers in 
parenthesis are the values of {X,A). 
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FIG. 10. The vertical velocity of the RG flow V2L — Vl against Vl for the N = 2 chiral model 
and the STA-XY model. "0(4)" denotes the 0(4)-symmetric fixed point. 
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FIG. 11. Histogram of energy per spin of the 1/3,2 niodel with the size L = 64 and L = 80 at 
K = 0.217685. 
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FIG. 12. RG flow of the = 3 chiral cf)^ model and the STA-Heisenberg model at the strong 
anisotropy region. The numbers in parenthesis are the values of (A, A) of the chiral 0^ model. 
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FIG. 13. RG flow of the N = 3 chiral cf)^ model at the weak anisotropy region. The numbers in 
parenthesis are the values of {\,A). 
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FIG. 14. The vertical velocity of the RG flow V2L — Vl against Vl for the N = 3 chiral model 
and the STA-Heisenberg model. "0(6)" denotes the 0(6)-symmetric fixed point. 
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FIG. 15. The RG flow of tlic N = 8 chiral cp'^ model. The numbers in parenthesis show the 
values of {X,A). CI and C2 denote the stable and unstable RG fixed point, respectively. 
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FIG. 16. The vertical velocity of the RG flow V2L — Vl against Vl near the stable chiral fixed 
point Ci of the N = 8 chiral cf)^ model. 



21 



A 


A 


K/N 


CiNT 


0.4 
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0.18 
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0.4 


0.8 
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0.4 


1.4 
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1.4 


0.4 
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2.1 


2.0 


0.0 


0.2880 


0.12 


2.0 


0.1 


0.3095 


0.12 


2.0 


0.2 


0.3292 


0.12 


2.0 
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0.3480 
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0.32 
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2.0 
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0.4510 


0.40 
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0.00 


0.2385 


0.14 


20.0 


0.02 


0.2415 


0.14 


20.0 


0.06 


0.2450 


0.21 


20.0 


0.09 


0.2460 


0.35 


STA 


0.7726 


0.10 



TABLE I. Parameters A, A, and K (divided by A^) at whicli simulations were carried out for 
the N = 2 chiral model and the STA-Heisenberg model. Physical quantities were observed for 
every Cint x L'^ MCS. 
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TABLE II. Parameters A, A, and K (divided by A'^) at which simulations were carried out for 
the N = 3 chiral (p^ model and the STA-Heisenberg model. Physical quantities were observed for 
every Cint x L"^ MCS. 
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TABLE III. Parameters A, A, and K (divided by N) at which simulations were carried out for 
the N = 8 chiral model. Physical quantities were observed for every Cjnt x L'^ MCS. 
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